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I .  INTRODUCTION 

Thermoacoustic  heat  transport  is  a  process  through  which 
an  acoustic  field  generates,  or  inversely  is  generated  from, 
a  flow  of  heat.  Like  every  engine,  a  thermoacoustic  engine  can 
be  configured  as  either  type  of  classical  heat  engine,  a  heat 
pump (or  refrigerator)  or  a  prime  mover.  The  purpose  of  this 
thesis  is  the  computer  implementation  of  a  general  formulation 
of  thermoacoustics  developed  by  Arnott  et  al  [Ref.  1]  using 
Matlab.  In  particular,  we  find  the  theoretical  values  of  the 
quality  factor  and  the  resonance  frequency  of  a  prime  mover  as 
a  function  of  the  temperature  difference  applied  across  the 
prime  mover  stack.  The  program  determines  the  impedance  and 
normalized  pressure  distribution  along  the  prime  mover 
beginning  from  a  known  value  of  specific  acoustic  impedance 
and  normalized  pressure  amplitude  at  one  end.  It  finds  the 
resonance  frequency  and  quality  factor  through  calculation  of 
the  pressure  amplitude  at  the  opposite  end  as  a  function  of 
frequency.  The  result  of  the  computer  implementation  are 
compared  with  measurements  made  with  open  and  closed  ended 
prime  mover  [Ref.  2,3],  and  with  the  results  of  a  standing 
wave  analysis  of  prime  movers  [Ref.  4]. 


II .   THEORY 

The  goal  of  this  chapter  is  to  describe  the  basic 
geometry  of  a  prime  mover  and  to  summarize  the  important 
points  of  Arnott's  method.  Parallel  with  this  we  explain  how 
this  theory  is  implemented  in  the  computer  program.  The  reader 
is  referred  to  [Ref.  1]  for  full  details  of  Arnott's  method. 

A.   ARNOTT'S  METHOD. 

The  necessary  background  for  this  thesis  can  be  explained 
with  the  help  of  Fig.l  which  shows  the  basic  construction  of 
the  thermoacoustic  prime  mover.  It  constists  of  a  circular 
cross  section  acoustic  resonator  ,  rigidly  capped  at  both 
ends.  Situated  within  the  resonator  are  the  thermoacoustic 
elements  consisting  of  an  ambient  heat  exchanger,  a  prime 
mover  stack  and  a  hot  heat  exchanger.  In  our  prime  mover  both 
the  stack  and  heat  exchangers  are  parallel  plates  spaced  by  a 
few  thermal  penetration  depths.  Complete  details  of  the  prime 
mover  construction  are  given  in  [Ref.  2,3,4].  The  traverse 
coordinates  of  the  prime  mover  are  taken  to  be  x  and  y,  and 
the  longitudinal  coordinate  is  z.  The  hot  and  cold  sections 
are  circular  ducts  open  at  the  heat  exchanger  end  and  open  or 
closed  at  the  other  end.  The  prime  mover  is  divided  into 
either  5  or  7  sections  depending  on  the  boundary  conditions 
[Fig.  1,2].  The  stack  is  divided  into  11   subsections. 
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Figure  1.  Rigid  end  Prime  Mover 
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Figure   2.    Open   end   Prime   Mover 


The  computer  program  finds  the  normalized  acoustic 
pressure  and  the  specific  acoustic  impedance  in  all  parts  of 
the  prime  mover  as  a  function  of  frequency.  The  pressure  as  a 
function  of  frequency  gives  the  frequency  response.  From  the 
frequency  response  we  can  determine  the  quality  factor  (Q)  and 
the  resonance  frequency  (f0)  .  To  do  that  we  assume  a  linear 
temperature  gradient  across  the  stack  via  the  hot  and  ambient 
(cold)  heat  exchangers.  The  stack  is  divided  into  eleven 
subsections.  A  constant  temperature  is  assumed  in  each 
subsection  of  the  stack.  In  all  others  parts  of  the  prime 
mover  we  assume  that  the  temperature  is  constant  and  equal  to 
that  of  the  nearest  heat  exchanger.  The  ambient  temperature 
(T0)  is  taken  to  be  a  function  of  z  (T0(z))  only  in  the  stack 
area . 

The  program  starts  with  specific  acoustic  impedance  at 
the  right  end.  There  are  two  boundary  conditions  of  interest: 
rigid  and  open.  The  specific  acoustic  impedance  of  a  rigid  end 
is : 


=  (l+i)^ 


2  (Y  -DPoCN 


P  0C 


v^ 


(1) 


as  described  in  [Ref.  5:p.529].  In  the  above  equation  Npr  is 
the  Prandtl  number,  cp  is  the  isobaric  heat  capacity  per  unit 
mass,  y  is  the  ratio  cp/cv where  cv  is  the  constant  volume  heat 


capacity  per  unit  mass,  T|  is  the  viscosity  and  p0  represents 
the  ambient  value  of  the  density.  The  specific  acoustic 
impedance  of  an  open  end  is: 


Z=cp  0Jca(^-0.6i),  (2) 

4 


where  a  is  the  radius  of  the  tube  and  k=0)/c  is  the  wavenumber, 
CO  is  the  angular  frequency  and  i=(-l);*  [Ref.  6:p.202]. 

The  acoustic  pressure  is  specified  at  the  same  end. 
Because  the  Q  is  independent  of  the  absolute  value  of  the 
acoustic  pressure  in  linear  theory,  the  specified  pressure  is 
only  a  relative  or  normalized  pressure.  Next  we  find  the 
pressure  and  the  impedance  at  the  entrance  to  the  hot  heat 
exchanger  using  the  Rayleigh' s  impedance  translation 
theorem  [Ref .1,5] .  This  is  one  difference  between  Arnott's  and 
other  methods.  That  is,  other  methods  find  the  fluid  variables 
by  integrating  the  governing  equations  along  the  prime  mover. 
The  impedance  translation  theorem  state  that:  within  any 
homogeneous  layer,  with  intrinsic  (or  characteristic)  specific 
acoustic  impedance  Zint,  the  local  specific  impedance  Z(z-l)  at 
z-1    is   related   to   that    at    z  by    [Ref.    5 :p. 139] : 


Z(y)  cos  (Jcl)  -iZlntsin(JcI) 
ZJntcos  (kl)  -iz(y)  sin(icl) 


z{z-i)=zint„J'  _  ,;,/  ,_/nvt  .   ;, ,;  .  o) 


Likewise,  the  acoustic  pressure  at  the  entrance  to  the  hot 
heat  exchanger  is  found  through  the  pressure  translation 
theorem: 

P1(z-l)=P1(z){cos(kl)-i[-^-]sin(kl)),  (4) 

Z(z) 

where  Z(z)  is  the  specific  impedance  at  z,  Pjfz)  is  the  first 

order  value  of  pressure,  i.e,  the  acoustic  pressure,  and  Zinc 

is: 


Z   -    Po"  (5) 

int     [C1F{K  )k]  ' 

In  the  above  equation,  Q  represents  the  porosity  and  is  equal 
to  NA,  where  A  is  the  cross  sectional  area  of  single  pore  and 
N  is  the  number  of  pores  per  unit  area,  k  is  the  complex  wave 
number  in  the  pore  and  given  by: 


(Jc(A  ,A  r))2  =  i^i   1_^(Y  -(Y  -1)F(A  T))  .        (6) 
C2  F(A  ) 


The  F  (X)     and  F  (XT)     factors  arise  from  the  equation  for 
the  z  component  of  the  acoustic  velocity: 


u,U,y,z)-F'*-^>diVZ)  .  (7) 

z  iup0    dz 


The  F(x,y;A.)  satisfies  the  following  differential  equation 
[Ref .  1] : 


R2 


F(x,y;X   )  +  (  -2-  )  V?.F(x,y;  A  )=1,  (8) 

ih  2 

where  the  Laplacian  operator  V2T  is  given  by: 

v°t=J_±+jLL.  (9) 

T   Bx2    By2 

In  our  case  where  the  pores  of  the  stack  are  parallel  plates, 
the  F  depends  only  on  y  and  the  shear  wave  number  X  or  the 
dimensionless  thermal  disturbance  number  X.T  which  are  given 
by: 


x  =R(pj^-)i  do) 


and 


A  r=A  Npr. 


The  F(y;A,)  is  given  by 


(11) 


cosh(/=T—  Z) 

F(y,k   )=1- 1T~-  (12) 

cosh(/=T— ) 
2 


In  the  above  equation  a  is  the  half  separation  distance 
between  the  parallel  plates  of  the  stack  or  heat  exchanger. 
The  pore  average  of  F(y;A.)  for  heat  exchangers  and  stack,  is 
given  by: 

F(A  )=l-(^-yTI)t&nh{yf=ri^-)  .  (13) 

A  2 


Finally  we  use  the  boundary  layer  approximation  in  the 
circular  regions  of  the  prime  mover: 

limx^F(A  )=l-(5  /R)  (l+i)  ,  (14) 

In  the  above  relation  R  is  twice  the  hydraulic  radius  of  the 
pore,  or  twice  the  ratio  of  the  transverse  pore  area  to  the 
pore  perimeter.  Actually  for  heat  exchangers  and  stack  with 
parallel  plates,  R  is  the  separation  distance  (2a)  between  the 
plates  and  for  the  circular  parts  it  is  the  radius  of  the 
tube.  8  is  the  viscous  or  thermal  penetration  depth. 

Once  Px  and  Z  are  known  at  the  entrance  to  heat  exchanger, 
we  need  the  values  just  inside  the  heat  exchanger.  Using  the 
pressure  and  the  impedance  translation  theorem  we  have  a 
solution  up  to  the  stack  where  a  linear  temperature  gradients 
exists.  In  the  heat  exchanger  three  things  are  different  from 
the  circular  parts.  These  are:  the  wavenumbers  in  the  heat 


exchangers  are  complex,  the  porosity  now  is  NA,  and  F  (X)    as  we 
mention  above  is  given  by  equation  (13) . 

The  point  now  is  how  to  find  pressure  and  impedance  in 
the  stack.  In  this  case  Arnott's  method  uses  the  following 
system  of  the  first  order  differential  equations  that  can  be 
solved  using  Range-Kutta  methods: 


dzJz)  =ik(z)Zint(z)  (i-  z(z)2  )+2a  (z)Z(z)  ,     (15) 


and 


dPi(z)   -^,     s*        ,  ,  Pl<z> 

— ij =ik  ( z)  Zint  ( z)     *    , 

dz  xnc         Z(z) 


where  a(z)  is 


Toz  is  the  temperature  gradient  given  by 


Oz 


dz 


(3  is  the  thermal  expansion  coefficient: 


10 


(16) 


f  (A  r) 

(A.xj-lSai  gfiLiL  ],  (17) 

2      1-Wp, 


T„-dT°(z),  (18) 


3  = 


e  r  p  (i9 

Po 


In  our  program,  this  system  of  equations  is  solved  with  the 
Runge-Kutta-Fehlberg  method  for  every  subsection  of  the  stack. 
In  so  doing,  the  impedance  and  the  pressure  distribution  in 
the  stack  is  obtained.  Now  using  the  pressure  and  the 
impedance  translation  theorem  we  find  the  pressure  and  the 
impedance  in  the  rest  of  the  prime  mover. 

Knowing  the  pressure  at  the  left  end  we  can  calculate  the 
quality  factor,  the  resonance  frequency  and  the  maximum 
pressure  amplitude.  This  is  done  by  minimizing  the  following 
equation  [Ref.  8]: 


<|p",'|fv 

— — *± -amplit2,  (20) 


f     to        O2 


where  Q  is  the  quality  factor,  |P(1) I  is  the  maximum 
amplitude  of  pressure  at  the  left  end,  f0  is  the  resonance 
frequency,  amplit  is  the  calculated  amplitude  of  pressure  at 
the  left  end  and  is  a  function  of  the  frequency  f.  In  the 
above  relation  there  are  two  known  variables  f  and  amplit  and 
three  unknowns  (Q,f0,P(l))  that  we  want  to  find.  To  find  these 
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unknowns,  the  3-parameter  least  squares  technique  is  used  with 
the  help  of  simplex  method. 

Up  to  this  point,  Q  and  f0  have  been  determined  only  by 
considering  the  pressure  amplitude  at  the  left  end.  We  have 
not  worried  about  matching  the  acoustic  impedance  at  the  left 
end.  Impedance  matching  allows  us  to  refine  the  values  of  Q 
and  f0  through  the  complex  eigenf requency .  Using  this 
technique,  one  assumes  the  angular  frequency  is  complex  and 
given  by: 


u=2nf0(l-i.  (21) 

Substituting  the  initial  value  of  f0  and  Q  determined  from  the 
frequency  response,  into  this  equation  gives  an  initial  guess 
at  the  complex  eigenf requency  for  the  system.  This  frequency 
is  used  to  calculate  the  impedance  throughout  the  prime  mover 
as  explained  earlier.  The  value  of  impedance  in  the  gas  at  the 
left  end  is  compared  to  the  value  for  the  boundary  condition 
(i.e.  the  impedance  of  a  rigid  end  with  thermal  loss).  The 
difference  between  the  computed  and  actual  impedance  is  used 
to  adjust  the  eigenf requency .  The  impedance  is  computed  with 
the  new  eigenf requency  until  both  real  and  imaginary  parts  of 
the  impedance  matche  at  the  left  end.  The  final  values  of  f0 
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and  Q  are  determined  from  the  real  and  imaginary  parts  of  the 
eigenf requency  that  best  matche  the  impedance. 

B.   IMPLEMENTATION  OF  ARNOTT' S  METHOD 

The  program  uses  functions  that  already  exist  in  Matlab 
where  possible.  For  the  Runge-Kutta-Fehlberg  5th  order  method 
we  modified  the  existing  function  in  Matlab  to  integrate  with 
initial  values  of  the  integration  variable  greater  than  the 
final.  This  is  necessary  because  we  put  the  origin  at  the  left 
end  (z  =  0)  and  we  begin  our  calculations  from  the  right  end 
(z=L)  .  One  reason  that  we  use  this  particular  Runge-Kutta 
method  is  its  accuracy.  The  Runge-Kutta-Fehlberg  computes  two 
Runge-Kutta  estimates  for  the  value  yn+1  but  of  different 
orders  of  error.  The  global  error  in  this  numerical  method  is 
0(h5)  and  the  local  error  is  0(h6)  [Ref.  7]. 

To  do  the  minimization,  we  use  the  function  fmins  with 
the  help  of  simplex  numerical  method.  The  fmin  function  with 
simplex  method  is  used  to  match  the  impedances  at  the  left  end 
using  the  technique  of  least  squares. 
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III.  THE  MATLAB  PROGRAM 

The  goal  of  this  chapter  is  to  provide  an  overview  of  the 
structure  of  the  Matlab  program  that  implements  Arnott's 
method.  Furthermore  it  describes  all  the  functions,  special 
tricks  and  limitations  of  this  code.  It  also  explains  how  to 
run  the  program. 

A.   OUTLINE  OF  THE  MATLAB  PROGRAM 

The  program  is  divided  into  three  parts.  In  the  first 
part  the  variables  are  specified  for  every  section  of  the 
prime  mover.  They  are:  the  type  of  section  (open,  heat 
exchanger,  stack) ,  the  number  of  subsections  that  exists 
inside  the  section,  the  temperature  at  the  left  and  right  ends 
of  the  section,  twice  the  hydraulic  radius  of  the  section,  and 
the  porosity  of  the  section.  Also,  specified  are  the  frequency 
range  and  the  number  of  intervals  in  this  range. 

The  second  part  is  the  main  program.  It  finds  the 
impedance  and  the  pressure  distribution  in  every  section  of 
the  prime  mover  for  each  frequency  within  the  assumed  range. 
To  accomplish  this  it  begins  at  the  right  end  of  the  tube 
where  it  calculates  the  impedance  of  the  end,  rigid  or  open, 
and  assigns  an  arbitrary  value  to  the  pressure.  Continuing,  it 
calculates  using  the  translation  theorems,  the  pressure  and 
the  impedance  in  all  other  sections  except  in  the  stack.  In 
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the  stack,  where  there  exists  a  temperature  gradient,  the 
program  calls  the  function  "trode45"  (a  modified  version  of 
the  function  of  Matlab  "ode45")  to  integrate  from  the  right 
end  to  the  left  end  of  the  stack.  This  function  solves 
differential  equations  using  the  5th  order  Runge-Kutta 
numerical  method  with  Fehlberg  coefficients. 

The  last  part  of  the  Matlab  program  uses  the  calculated 
pressure  at  the  left  end  from  the  previous  part  to  find  the 
maximum  amplitude,  the  resonance  frequency  and  the  quality 
factor.  These  parameters  are  found  from  a  least  square  fit  to 
equation  (20) .  The  minimization  is  performed  with  the  matlab 
function  "fmins"  with  "options  1,2,3,14".  Finally,  the  Q  and 
f0  estimates  are  used  as  a  starting  points  for  the  final 
calculation  of  the  quality  factor  and  the  resonance  frequency 
using  the  complex  eigenf requency  method.  This  refinement 
usually  is  a  small  correction  (in  the  fourth  decimal  place)  to 
initial  guesses.  Finally,  the  program  plots  1/Q  and  f0  versus 
temperature  difference  and  saves  the  results  in  a  data  file. 


B.   LIMITATIONS  OF  THE  PROGRAM,  SPECIAL  POINTS. 

There  are  two  programs,  one  for  the  open  end  prime  mover 
and  one  for  the  rigid  end  tube.  To  run  the  open  end  progam, 
one  needs  to  run  "arffl.m"  with  the  functions  "error, m- 
erro3 .m-argo3f .m" .  To  run  the  rigid  end  program  one  needs  to 
run  "argo2.m"  with  the  functions  "error .m-erro3 .m-argo3 .m" . 
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The  difference  in  the  two  progams  is  that  one  uses  "argo3.m" 
and  the  other  uses  "argof 3 .m" .  Before  running  one  program  or 
the  other,  one  needs  to  switch  these  functions  in  the 
"erro3.m"  function. 

Futhermore,  there  are  some  factors  that  influence  the 
operation  of  the  program.  One  important  step  is  the  choice  of 
the  initial  range  of  the  frequencies.  To  begin  the  program  one 
needs  to  know  the  approximate  resonance  frequency  as  well  as 
which  longitudinal  acoustic  mode  is  to  be  used.  Initially  the 
frequency  range  is  f0  ±  20Hz.  After  the  program  runs  for  the 
first  value  of  AT  the  frequency  range  is  automatically 
adjusted  to  be: 

Af=-i.  (22) 

151 

Another  important  factor  is  of  course  the  number  of 
points  between  the  maximum  and  minimum  frequency  that  the 
program  uses.  Usually  this  number  is  500  but  sometimes  when 
the  values  of  quality  factor  gets  large  this  may  have  to  be 
doubled  (1000  points)  .  In  these  cases,  it  may  also  be 
necessary  to  use  steps  of  5  degrees  Kelvin  instead  of  the 
usual  10  degrees  Kelvin. 

Finally,  another  important  point  concerns  the  "counters" 
that  are  used  by  the  program.  The  most  important  counter  is 
the  "flag".  If  somebody  wants  to  run  the  program  for  a  case 
where  the  quality  factor  increases  with  AT,  one  needs  to  put 
the  initial  value  of  the  "flag"  counter  equal  to  zero. 
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Otherwise,  it  should  be  set  equal  to  one.  The  purpose  of  this 
counter  is  to  force  the  program  to  scan  for  negative  values  of 
the  inverse  quality  factor  at  the  proper  time. 
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IV.  VALIDATION  AND  RESULTS 

In  this  chapter  we  compare  the  results  of  the  program  to 
experimental  results  for  two  different  prime  movers  [Ref.  2 
and  3]  as  well  as  with  the  results  of  the  standing  wave 
analysis  by  Atchley  [Ref  4]. 

The  first  results  are  for  a  rigid  end  prime  mover 
described  in  Ref.  2  and  4.  Comparison  of  both  1/Q  and 
resonance  frequency  are  made.  Figures  3  through  5  show  1/Q  vs 
AT  (in  degrees  K)  for  three  different  mean  gas  pressures.  The 
gas  is  helium.  In  these  figures,  the  open  circles  represent 
measured  values,  the  *'s  represent  the  results  of  the  standing 
wave  analysis  and  the  solid  line  represents  the  result  of  the 
Matlab  program.  The  overall  agreement  between  the  present 
analysis  and  the  measurements  is  good,  but  not  as  good  as  for 
the  standing  wave  analysis.  The  present  analysis  agrees  best 
at  low  values  of  AT.  The  reasons  for  the  increased  discrepancy 
at  higher  values  of  AT  are  not  understood.  One  would  expect 
the  present  analysis  to  agree  better  with  measurements  than 
the  standing  wave  analysis.  This  is  expected  because  the 
present  analysis  makes  fewer  assumptions.  In  the  particular, 
it  does  a  better  job  of  taking  into  account  changes  in  the 
acoustic  pressure  amplitude  in  differents  parts  of  the  prime 
mover.  The  standing  wave  analysis  ignores  any  changes  in 
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pressure  amplitude  due  to  changes  in  cross  sectional  area  of 
the  prime  mover. 

Figures  6  through  8  show  results  of  resonance  frequency 
(in  Hz)  vs  AT  for  the  same  three  mean  pressures.  The  agreement 
between  the  present  analysis  and  the  measured  frequencies  is 
within  1%  in  all  cases.  The  standing  wave  analysis  also  agrees 
with  the  measured  values  to  about  the  same  accuracy.  However, 
the  two  methods  show  different  depedences  on  AT.  This  is 
because  1)  the  standing  wave  analysis  considers  only  the 
temperature  dependence  of  the  sound  speed,  whereas  the  present 
analysis  includes  the  effects  of  dispersion  due  to  the 
boundary  layers,  and  2)  the  standing  wave  analysis  ignores  the 
finite  impedance  of  the  rigid  ends.  Of  these  two  differences 
the  first  is  the  more  important.  This  results  in  different 
curvatures  of  the  two  theoretical  lines. 

The  second  set  of  comparisons  is  made  for  an  open  ended 
prime  mover  used  in  Che's  thesis  research  [Ref.  3].  The 
unusual  property  of  this  prime  mover  is  that  while  1/Q 
decreases  with  AT  for  the  first  and  third  modes  (as  is  usual 
for  a  prime  mover) ,  1/Q  increases  with  AT  for  the  second  mode. 
This  behavior  provides  a  good  test  for  theoretical  models.  The 
results  for  1/Q  vs  AT  are  shown  in  figures  9  through  11  for 
the  first,  second  and  third  modes,  respectively.  It  is  seen 
that  the  present  analysis  and  the  standing  wave  analysis  are 
in  close  agreement  in  the  first  and  second  modes.  The 
agreement  diminishes  slightly  at  the  third  mode.  The  agreement 
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with  the  measured  values  is  worse  than  before.  However,  this 
is  most  likely  due  to  contamination  of  the  helium  with  air  and 
water  vapor  as  discuseed  in  Che's  thesis.  Another  interesting 
aspect  of  the  analysis  is  the  prediction  of  an  extremum  in  1/Q 
for  the  second  and  third  modes. 
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Figure  3.  Graph  of  1/Q  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  170  kPa  mean  pressure. 
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Figure  4.  Graph  of  1/Q  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  376  kPa  mean  pressure. 
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Figure  5.  Graph  of  1/Q  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  500  kPa  mean  pressure. 
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Figure  6.  Graph  of  f0  (Hz)  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  170  kPa  mean  pressure. 
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Figure  7.  Graph  of  f0  (Hz)  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  376  kPa  mean  pressure. 
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Figure  8.  Graph  of  f0  (Hz)  vs  AT  (K)  for  the  closed  end  prime 
mover  filled  with  helium  at  500  kPa  mean  pressure. 
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Figure  9.  Graph  1/Q  vs  AT  (K)  for  the  first  mode  of  the  open 
ended  prime  mover  filled  with  helium  gas  at  101  kPa. 
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Figure  10.  Graph  of  1/Q  vs  AT  (K)  for  the  second  mode  of  the 
open  ended  prime  mover  filled  with  helium  gas  at  101  kPa. 
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Figure  11.  Graph  of  1/Q  vs  AT  (K)  for  the  third  mode  of  the 
open  ended  prime  mover  filled  with  helium  gas  at  101  kPa 
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Figure  12.  Graph  of  f0  (Hz)  vs  AT  (K)  for  the  first  mode  of 
the  open  ended  prime  mover  filled  with  helium  gas 
at  101  kPa. 


30 


540 


520 


o 

c 


500 


S"480 


o 

§460 

c 
o 

CO 

to 

^  440 


420 


400 


1 

■1 

1                 1                 1 

Ambient  pressure  101  kF 

a 

i 

Seconc 

mode 

20  40  60  80         100        120 

Negative  DeltaT 


140    160 


180 


Figure  13.  Graph  of  f0  (Hz)  vs  AT  (K)  for  second  mode  of  the 
open  ended  prime  mover  filled  with  helium  gas  at  101  kPa. 


31 


900 

880 

860 
o 

SJ840 
cf 

LL 

0,820 

o 

c 
CO 

§800 

0) 

cc 
780 

760 

740 

i 

i               i               i 

Amhiont  nrfft^iirn  1f)T  kPn       ... 

"v.                Third  rriode 

^^_i 

i 

>               > 

)              20             40             60             80            100           120           140 

Negative  DeltaT 
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the  open  ended  prime  mover  filled  with  helium  gas 
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V.  SUMMARY 

The  results  of  this  thesis  can  be  summarized  as  follows. 
Arnott's  method  was  implemented  in  a  Matlab  progam.  The 
impedance  and  the  pressure  distribution  were  calculated  along 
the  prime  mover.  From  the  frequency  response  at  the  left  end 
and  complex  eigenf requency,  the  quality  factor  and  the 
resonance  frequency  was  estimated  for  both  cases  rigid  and 
open  end  prime  movers.  The  results  of  the  analysis  were 
compared  with  those  of  a  standing  wave  analysis  and  with 
experimental  results.  In  the  most  of  the  cases  the  results 
agree  well.  A  future  investigation  can  be  the  examination  of 
the  output  work  of  the  prime  mover. 
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APPENDIX  A.  THE  MATLAB  PROGRAM  FOR  THE  OPEN  END  PRIME  MOVER 


%  Program  for  prime  mover 

%  PART  1   (CHANGE) 
clear; 
clg; 
coun=0; 
pt=0; 

qual  =  zeros  (45, 1) ; 
inqual  =  zeros  (45, 1) ; 
frO  =  zeros  (45, 1) ; 
deltaT=zeros  (45, 1) ; 
point=0; 
flag=l; 

slope=sum (' negative' ) ; 
load  open.dat 


%  counter 

%  counter 

%  quality  factor 

%  inverse  quality  factor 

%  resonant  frequency 

%  temperature  difference 

%  indicator 

%  indicator 

%  standing  wave  analysis  data 


Main  Loop  for  evrey  temperature  difference  (itrgh) 
for  itrgh=293:-5:158 
coun=l+coun 

numfre=1000;  %  number  of  frequencies 

s   establish  some  often  used  constants 

s=5  ;  %  number  of  sections 

j=sqrt (-1) ; 
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npr=2/3* (1+eps) ; 

gamma=5/3* (1+eps)  ; 

cp=(2.5*8.314  3/ (4.e-3) )* (1+eps) ; 

ksolid=. 16* (1+eps) ; 

sqri=( (l+j)/sqrt (2) ) * (1+eps) ; 

pl=ones (18, numf re) . * (10^ (-8) ) ;  %  presure  inside  the  tube 

f req=zeros (1,  numf re) ;    %  frequncies 


w=zeros (1, numf re) ; 
teren=zeros (1,  s)  ; 
numblay=zeros (1,  s)  ; 
lengsec=zeros (1,  s)  ; 
tempR=zeros (1,  s)  ; 
tempL=zeros ( 1 ,  s )  ; 
poros  =  zeros (1,  s)  ; 
ratio=zeros (1,  s)  / 
%  SET  VALUES 

f remin=864; 

f remax=904; 
termin=sum (abs (' free' ) ) ; 
ampres= (1 . 01e+5) * (1+eps)  ; 

dramp= (1 . e-16) ; 


%  angular  frequencies 

%  end  of  sections 

%  layers  of  sections 

%  lengths  of  sections 

%  rigth  temperature  of  sections 

%  left  temperature  of  sections 

%  porocity  of  stack 

%  ratio 


%  minimum  frequency 

%  maximum  frequency 

%  end  of  the  last  section 

%  ambient  presure 

%  driver  pressure 


%  section  "1" 

terenl  =  ' opentu'  ; 

teren (1,1) =sum (abs (terenl)  )  ; 

numblay (1, 1) =1; 
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lengsec 
tempR(l 
tempL (1 
ratio  (1 
poros  (1 
ares  (1, 

;  section 
teren2= 
teren (1 
numblay 
lengsec 
ternpRd 
tempL (1 
ratio  (1 
poros (1 

%  section 
teren3= 
teren (1 
numblay 
lengsec 
ternpRd 
tempL (1 
ratio (1 
poros (1 

%  section 
teren4= 


1, l)=(66.25e-2) * (1+eps) ; 

1)=293*  (1+eps) ; 

1)=293*  (1+eps)  ; 

l)  =  (1.917e-2)  *  d+eps)  ; 

1)=1* (1+eps) ; 

) =pi*ratio (1, : ) * (1+eps) ; 

2" 

hexch' ; 

2) =sum  (abs (teren2) ) ; 

1,2)=1; 

l,2)=(.82e-2) * (1+eps) ; 

2)=293*  (1+eps) ; 

2) =293*  (1+eps) ; 

2)=(5.08e-4) * (1+eps) ; 

2)=. 666*  (1+eps) ; 

"3" 

stack' ; 

3) =sum (abs (teren3) ) ; 

1,3)=11; 

l,3)=(2.59e-2) * (1+eps) ; 

3) =itrgh* (1+eps) ; 

3) =293*  (1+eps) ; 

3)  =  (8.636e-4) *  (1+eps) ; 

3)=. 89473* (1+eps) ; 
ii  ^  ii 

hexch' ; 
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teren  (1,4) =sum (abs (teren4) ) ; 
numblay (1, 4)=1; 

lengsec  (1,4)  =  (. 82e-2) * (1+eps) ; 
tempRCL,  4)  =itrgh*  (1+eps)  ; 
tempL (1, 4) =itrgh* (1+eps) ; 
ratio(l,4)=(5.08e-4) * (1+eps) ; 
poros  (1, 4) =. 666* (1+eps) ; 
%  section  "5" 

teren5=' opentu' ; 

teren  (1,5) =sum (abs (teren5) ) ; 

numblay (1, 5) =1; 

lengsec (1,5)  =  ( . 6852) * (1+eps)  ; 

tempR(l, 5) =itrgh* (1+eps) ; 

tempL  (1, 5) =itrgh* (1+eps) ; 

ratio(l,5)=(1.917e-2) * (1+eps) ; 

poros (1, 5) =1* (1+eps) ; 

%  PART  2 
%   MAIN  PROGRAM 
%  FREQUENCY  DISTRIBUTION 

numtot=sum (numblay) +1; 

pt=pt+l; 

if  pt~=l 

df=frO  (coun-1) /abs (qual (coun-1)  )  ; 

fremin=frO (coun-1) -df ; 
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fremax=frO (coun-1) +df ; 
end 

i=l :numf re; 
if  i~=l 

f req=f remin; 
else 

f req (1,  i) =f remin+ (fremax-f remin)  . * (i . / (numf re-1) ) ; 
end 
w=2*pi*freq; 
%  GET  THE  SPECIFIC  ACOUSTIC  IMPENDANCE  AND  PRESSURE  AT  ALL 

POINTS. 
%  START  AT  THE  RIGHT  AN  MOVE  AT  THE  LEFT, 
dens (1, numtot) =ampres*4. Oe-3  / (tempR (1 , s) *8 . 3143) ;  %  density 
visc(l,numtot) =1.887e-5* (tempR (1, s) /273.15) A (.65  67) ;% 

viscosity 
kgas (1, numtot ) =visc (1, numtot ) *cp/npr;  %  termal  conductivity 
sspeedd,  numtot)  =972.  8*sqrt (tempR(l,s) /273.15) ;  %   speed 

%   isothermal 
typel  =  /  free'  ; 
type2  =  / rigid'  ; 

if  termin==sum (abs (type2) ) 
%   IMPEDANCE  OF  RIGID  TERMINATION. 

fac (numtot, 1 : numf re) =sqrt ( (dens (1, numtot ) . . . 

*s speed (1, numtot) ~2) ./(w.*visc(l, numtot) ) ) ; 
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z (numtot, 1 :numfre) = (1+ j) . * (dens (1, numtot) . *s speed (1, numtot ) . 
. . . *f ac (numtot, 1 :numf re) *sqrt (npr) ) / (sqrt (2) * ( gamma- 1) ) ; 
%  IMPEDANCE  OF  FREE  TERMINATION 
else 
z (numtot, 1 :numf re) =f ree (dens (1, numtot) , vise (1, numtot) , w, rati 
0(1,5)..., s speed (1, numtot) , gamma, npr) ; 
end 
%  IMPEDANCE  TRANSLATE  FOR  THE  OPEN  TUBE  OR  HEAT  EXCHANGER. 
fault=' stack' ; 
truel=' opentu' ; 
true2=' hexch' ; 
for  k=s:-l:l 
jup=0; 

jup=numtot-sum ( jup+numblay ( 1 ,  k  :  s )  )  ; 
if  teren (1, k) ~=sum (abs (fault) ) 
dens (1, jup)=ampres*4.0e-3  / (tempR (1, k) . *8 . 3143) ;  %  density 
vised,  jup)=1.8  87e-5.*  (tempR(l,k)  /273.15)  .  A  (.6567)  ; 

%  viscosity 
kgas  (1, jup) =visc (1, jup)  . *cp/npr;   %  termal  conductivity 
sspeedd,  jup)  =972  .  8  .  *sqrt  (tempR(l,k)  /273.15)  ;    %   speed 

%   isothermal 
%      GET  LAMBDA  AND  LAMBDA (T) 
lambda (jup, 1 :numf re) =ratio (1, k) . *sqrt (dens (1, jup ).*w. /vise (1 
/  jup) )  ; 
lambdt (jup, 1 : numf re) =sqrt (npr) . *lambda (jup, 1 :numf re) ; 
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%      GET  F(L),  F(L(T))  AND  WAVENUMBERS  FOR  OPEN  TUBE 

PARTS 
if  teren (1, k) ==sum (abs (truel) ) 
flam (jup, 1 :numf re) =1- (1+ j) . *sqrt (2) . /lambda ( jup, 1 mumfre) ; 

%  f(D 
f  lamt  (jup,  1  mumfre)  =1-  (1  + j)  .  *sqrt  (2)  .  /lambdt  (jup,  1  :numfre)  ; 

%  f  (Kt)) 
facl=  (1+  (gamma-1)  /sqrt  (npr)  )  /sqrt  (2)  ; 
ka  ( jup,  1  mumfre)  =  (w.  /s speed  (1,  jup)  )  .  .  . 

. * (1+  (1+j)  . *facl . /lambda (jup, 1 mumfre) ) ; 
else 
%      GET  F(L),  F(L(T))  AND  WAVENUMBERS  FOR  HEAT 
EXCHANGERS  TYPE  "SLIT". 

sqrmi= (1- j) /sqrt (2) ; 

ar  (jup,  1  mumfre)  =sqrmi  .  *lambda  (jup,  1  mumf  re)  ; 
argum  (jup,  1  mumfre)  =exp  (-ar  (jup,  1  mumf  re)  )  ; 
ar  ( jup,  1  mumfre)  =ar  ( jup,  1  mumf  re)  /2; 
ctanh ( jup, 1 : numf re) = (1 -argum (jup, 1 mumf re) ) . / (1+argum (jup, 1 : 
numfre) ) ; 
flam (jup, 1 mumfre) =1 -ctanh (jup, 1 mumfre) . /ar ( jup, 1 : numf re) ; 

%  f(D 
ar ( jup, 1 : numfre) =sqrmi*lambdt ( jup, 1 : numfre) ; 
argum (jup, 1 : numf re) =exp (-ar ( jup, 1 : numf re) ) ; 
ar ( jup, 1 : numf re) =ar ( jup, 1 : numfre) 12; 
ctanh ( jup, 1 : numfre) = (1-argum ( jup, 1 : numfre) ) . / (1+argum ( jup, 1 : 
numfre) ) ; 

40 


f lamt ( jup, 1 :numf re) =l-ctanh ( jup, 1 :numf re)  . /ar ( jup,  1 :numf re) ; 

%  f (1 (t)  ) 
low=w/sspeed (1, jup)  ; 
ka ( jup, 1 :numf re) =low. *sqrt ( (gamma- ( gamma- 1) . . . 

. *f lamt (jup, 1 :numf re) ) . /f lam ( jup, 1 :numf re) ) ; 
end 
%  TRANSLATION  THEOREM 
comden (jup, 1 :numf re) =dens (1, jup) . /flam (jup, 1 :numf re) ; 
zint (jup, 1 : numf re) =comden ( jup,  1 : numf re)  . *w. / (ka ( jup, 1 :numf re 
) . *poros (1, k) ) ; 

dsub  (1, k) =lengsec (1, k)  . /numblay (1, k) ; 
ar ( jup, 1 : numf re) =ka (jup, 1 : numf re) . *dsub (1, k) ; 
sn  ( jup, 1 : numf re) =sin(ar(jup,l: numf re) ) ; 
cs  ( jup, 1 : numf re) =cos (ar ( jup, 1 : numf re) ) ; 
ct (jup, 1 : numf re) =cs ( jup, 1 : numf re) . /sn ( jup, 1 : numf re) ; 
fac3 (jup, 1 : numf re) =zint (jup, 1 : numf re) . /z ( jup+1, 1 : numf re) ; 

z (jup, 1 : numf re) =zint (jup, 1 : numf re) . * (ct ( jup, 1 : numf re) . . . 
- j . *f ac3 (jup, 1 : numf re) )  . /  (f ac3 ( jup, 1 : numf re)  . *ct ( jup, 1 :numf r 
e)-j); 

pi (jup, 1 : numf re) =pl (jup+1, 1 : numf re) . * (cs ( jup, 1 : numf re) . . . 
-j.*fac3(jup,l : numf re) . *sn ( jup, 1 : numf re) ) ; 
else 
%  STACK  WITH  LINEAR  DISTRIBUTION  OF  TEMPERATURE  BETWEEN 
THE  ENDS  OF  THE  STACK 

toz=  (ternpRd,  k)  -tempL  (1,  k)  )  /lengsec  (1,  k)  ;  %  approximation 

dz/dt 
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ns=0; 

for  la=k+numblay (1, k) -1 : -1 : k 
ns=ns+l; 

tens=tempR(l, k) - (tempR(l, k) -tempL (1, k) ) *ns/numblay (1, k) ; 
tnsml=tempR(l,k)  -  (tempR(l,k)  -tempL(l,  k)  )  *  (ns-1)  /numblay  (l,k) 

tave (1, la) = (tens+tnsml) 12; 

dsub (1, k) =lengsec (1, k) /numblay (1, k) ; 

dens (1, la) =ampres*4 . Oe-3  / (tave (1, la) *8 . 3143) ;     %  density 

vise  (1,  la)=1.887e-5*  (taved,  la)  /273.  15)  A  (.6567)  ;  %  viscosity 

kgas  (1, la) =visc (1, jup) *cp/npr;  %  termal   conductivity 

sspeedd,  la) =972 . 8*sqrt  (tave  (1,  la)  /273.15)  ;  %   speed 

%   isothermal 

%      GET  LAMBDA  AND  LAMBDA (T) 
lambda (la, 1 : numf re) =ratio (1, k) . *sqrt (dens (1,1a) .*w./visc(l,l 
a)  )  ; 
lambdt (la, 1 : numf re) =sqrt (npr) *lambda (la, 1 : numf re) ; 

ar (la,  1 : numf re) =sqrmi* lambda (la, 1 : numf re) ; 

argum (la, 1 :numf re) =exp (-ar (la, 1 :numf re) ) ; 

ar (la, 1 : numf re) =ar (la, 1 : numf re) 12; 
ctanh (la, 1 : numf re) = (1 -argum (la, 1 : numf re) ) . / (1+argum (la, 1 :num 
fre)); 
flam (la,  1 : numf re) =1 -ctanh (la, 1 : numf re)  . /ar (la, 1 : numf re) ; 

%  f(D 
ar  (la,  1 :  numf  re)  =sqrmi*  lambdt  (la,  1 :  numf  re)  ; 
argum (la, 1 : numf re) =exp (-ar (la, 1 : numf re) ) ; 
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ar (la,  1 : numf re) =ar(la, 1 :numf re) /2; 

ctanh (la, 1 :numf re) = (1-argum (la, 1 : numf re) ) . / (1+argum (la, 1 :num 

fre) )  ; 

f lamt (la, 1 : numf re) =1 -ctanh (la, 1 : numf re)  . /ar  (la, 1 : numf re) ; 

%  f  (Kt)  ) 
low=w/sspeed (1, la) ; 
ka (la,  1 : numf re) =low. *sqrt ( (gamma- ( gamma- 1) . . . 

. *f lamt (la, 1 : numf re) ) ./flam (la, 1: numf re) ) ; 
zint  (la, 1 : numf re) =dens (l,la).*w... 
. /  (poros (1, k)  .*flam(la,l : numf re)  . *ka (la, 1 : numf re) ) ; 
alprim (la, 1 : numf re) =toz* (f lamt (la, 1 : numf re) . . . 
. /flam(la, 1 rnumfre) -1) . / (2*tave (1, la) * (1-npr) ) ; 


kal=ka (la, 1 rnumfre! 


i  . 


zintl=zint (la, 1 : numf re) . ' ; 

alprl=alprim (la, 1 : numf re)  . '  ; 

zetaO= (lengsec (1, k) - (ns-1) *lengsec (1, k) /numblay (1,  k) ) ; 

zetaf in= (lengsec (1, k) - (ns) * lengsec (1, k) /numblay (1, k) ) ; 

zO=z(la+l,l: numf re) . ' ; 

tol=l .e-4; 

[zeta, zout ] =t rode 4 5 (' tube' , zetaO, zetafin,  zO,tol,  kal,  zintl,  al 

prl); 

n=length (zeta) ; 

z  (la,  1 :  numf  re)  =zout  (n,  :  )  ; 

zdumy=zO; 

plO=pl (la+1, 1 :numfre) . ' ; 

tol=l.e-4; 
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[zeta,  pi  out  ]  =t  rode  4  5  ('  tubel'  ,zetaO,zetafin,plO,tol,kal,zintl 

,  zdumy) ; 

nl=length (zeta) ; 

pi (la, 1 :numfre) =plout (nl, : ) ; 

end 

end 

end 

%  MINIMIZATION  TO  GET  RESONANCE  FREQUENCY  AND  QUALITY 
FACTOR. 

tfun=- j*w. *z  (1, 1 :numf re) *dramp. /pi (1, 1 :numf re) ; 
for  a=l:s 

pi (a, 1 : numfre) =pl (a, 1 :numf re) . *tfun; 
end 

amplit=zeros  (1, numfre) ; 

amplit (1,1 : numfre) =abs (pi (1,1 : numfre) ) ; 
[maxamp, i ] =max (amplit) ; 
resf re=f req (i) 
amphaf =maxamp/sqrt (2) ; 
q=0; 

frehaf=0; 
if  point==0 

for  s=2: numfre 

if  amplit (s) >=amphaf   &  amplit (s-1) <=amphaf 
frehaf= (freq(s) +freq(s-l) ) 12; 
q=0 . 5*resf re/ (resf re-f rehaf ) 
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end 
end 
else 

for  s=2:numfre 

if  amplit (s) <=amphaf   &  amplit (s-1) >=amphaf 
ama=l 
frehaf= (freq(s) +freq(s-l) ) /2 
q=0 . 5*resf re/ (resf re-f rehaf ) 
end 
end 
end 
if  q~=0 

xO=zeros (3,1); 
xO (1, 1) =maxamp; 
xO  (2, 1) =resfre; 
x0(3,l)=q; 
options (1) =0; 
options (2) =1 .  e-4; 
options (3) =1 .e-6; 
options (14) =2000; 

est=fmins (' erro' , xO, options,  [ ] , amplit, f req) ; 
inqual (coun, 1) =est  (3, 1) A  (-1) ; 
frO (coun, 1) =est  (2, 1) ; 
deltaT (coun, 1) =abs (itrgh-293) ; 
qual (coun, 1) =est  (3, 1) ; 
end 
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if  coun~=l 

dif= inqual  (coun-1, 1) -inqual (coun,  1) ; 
if  dif<0  &  coun==2 

slope=sum (abs (' positive'  )  )  ; 
elseif  dif<0  &  coun==3 

slope=sum (abs ('positive'  )  )  ; 
end 

if  flag==0  &   slope~=sum (abs ('positive' ) ) 
if   dif<=0  |  q==0 
point=l; 
flag=l; 
for  s=2:numfre 

if  amplit (s) <=amphaf   &  amplit (s-1) >=amphaf 
frehaf=(freq(s) +freq(s-l)  )  12; 
q=0 . 5*resf re/ (resf re-f rehaf ) 
end 
end 
xO (1, 1) =maxamp; 
xO  (2, 1) =resf re; 
x0(3,l)=q; 
options (1) =0; 
options (2) =1 .e-4; 
options (3) =1 . e-6; 
options (14)=2000; 

est=fmins  (' erro' , xO, options,  [ ] , amplit, f req) ; 
inqual (coun, 1) =est  (3, 1) * (-1) ; 
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frO (coun, l)=est  (2, 1) ; 

deltaT (coun, 1) =abs  (itrgh-293) ; 

qual (coun, 1) =est  (3, 1)  ; 

end 
end 
end 

%  PART  THREE 
%  MINIMIZATION  USING  COMPLEX  ANGULAR  FREQUENCIES, 
options  (1) =0; 
options (2) =1 . e-4; 
options (3) =1 .e-4; 
options (14) =700; 

wl=2*pi*fr0  (coun, 1) * (1-j/ (2*qual (coun, 1)  )  )  ; 
y2=wl+wl*10.e-6; 
yl=wl-wl*10 . e-6; 

ella=fmin  (' erro3' , yl, y2, options, itrgh) ; 
frO (coun, 1) =real (ella) / (2*pi)  ; 

qual (coun, 1)=- real (ella)/ (2*imag (ella) ) ; 
inqual (coun, 1) =qual (coun, 1) A (-1)  ; 
end 
%     PLOTS  AND  STORE  THE  DATA   (CHANGE) . 

y=[ (inqual (1 :coun, 1) ) ' ;  (deltaT (1 :coun, 1) ) ' ; 
(frO  (l:coun, 1) ) ' ]; 

f id=fopen (' inqf3.txt'  ,  '  w+'  )  ; 
fprintf (fid, ' Inqual=%-9 . 6f   T=%-4.1f   RF=%-5 . 2f \n' , y) ; 
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load  compf2.dat 

figure 
plot (deltaT (1 :coun) , inqual (1 : coun) , compf2 ( : , 1)  , compf2 ( : , 2) , ' 
o' .  .  . 

,open(: , 1) , open ( : , 4) , '*' ) 

title (' Inverse  Quality  Factor  V.S  Temperature 
Difference . Open  End' ) ; 

grid 

xlabel ('Negative  Delta  "T"  '); 

ylabel ('1/Q' ) ; 

gtext ('Ambient  pressure  101  Kpa' ) 

gtext (' Experimental  Result  "o"') 
gtext  ('Computer  Result  "*•") 

gtext ('Third  Mode') 

print  -deps  inquf3 

figure 

plot (deltaT (1 :coun) , frO  (1 :coun) ) 

xlabel ('Negative  Delta  "T"  '); 

ylabel (' Resonance  Frequence'); 

title (' Resonance  Frequency  V.S  Temperature 
Difference .Open  End'); 

gtext ('Ambient  pressure  101  Kpa') 

gtext ('Third  Mode') 

grid 

print  -deps  frf3 
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%  FUNCTION  TRODE4  5.M 

function  .  .  . 

[tout, yout ] =t rode 4 5 (FunFcn, tO,tfinal,yO,tol,kal, zintl, alprl 

alpha=[l/4  3/8  12/13  1  1/2]'; 

beta=[  [1      0      0     0  0  0]/4 

[3       9       0     0  0  0]/32 

[  1932   -7200    7296     0  0  01/2197 

[  8341  -32832   29440  -845  0  01/4104 

[-6080   41040  -28352  9295  -5643  01/20520]'; 

gamma=[  [902880  0  3953664  3855735  -1371249  2770201/7618050 

[  -2090  0    22528  21970  -15048  -27360 ] /752400 


r 


%  initialization 

t=t0; 

hmax= (t-tf inal)  /5; 

hmin= (t-tf inal) /2000; 
h=(t-tfinal) /100; 

y=y0 ( : ) ; 

f=y*zeros (1,6); 

tout=t; 

yout=y . ' ; 

tau=tol*max (norm (y, ' inf '  )  ,  1)  ; 
%  Main  loop 

while (t>tf inal) & (h>=hmin) 
if  t-h>tfinal  ,h=t-tfinal;  end 
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%  Compute  the  slopes 
temp=feval (FunFcn, kal, zintl, alprl, t,  y) ; 
f ( : , 1) =temp( : ) ; 
for    j=l:5 
temp=feval (FunFcn, kal, zintl, alprl, t-alpha ( j) *h, y-h*f *beta ( : , 

j)); 

f  (:,  j  +  l)=temp(:) ; 
end 
%  Estimate  the  error  and  the  acceptaple  error 
delta=norm(h*f *gamma ( : , 2) , '  inf '  ) ; 
tau=tol*max  (norm  (y, '  inf  ) ,  1) ; 
%  Update  the  solution  only  if  the  error  is  acceptaple 
if  delta<=tau 
t=t-h; 

y=y-h*f *gamma ( : , 1 ) ; 
tout= [tout; t ]  ; 
yout= [ yout ; y . ' ] ; 
end 
%  Update  the  step  size 
pow=l/5; 
if  delta~=0.0 

h=min (hmax, 0 . 8*h* (tau/delta) Apow) ; 
end 
end 
if  (tfinaKt) 

disp(  'singularity  likely.') 
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t 

end 
%  FUNCTION  "TUBE.M"  GIVES  THE  DERIVATIVE  OF  IMPEDANCE 
function  dzdz=tube (kal, zintl, alprl,  zeta,  z) 
dzdz=( j. *kal.* zintl.* (1- (z. /zintl) . "2) +  2 . *alprl . *z) ; 

%  FUNCTION  "TUBE1.M"  GIVES  THE  DERIVATIVE  OF  PRESSURE 
function  dpdz=tubel (kal, zintl, zdumy, zeta,pl) 
dpdz= j . *pl . *kal . *zintl . /zdumy; 

%  FUNCTION  "ERRO.M"  FOR  "FMINS".  LEAST  SQUARS  TECHNIQUE. 

function  rmserr=erro (xO,  amplit,  f req) 
mamp=xO (1,1); 
fO=xO  (2, 1)  ; 
z=xO  (3, 1)  ; 
focl= (mamp*f 0) . / (z*f req) ; 

err=(  (foci.  "2)  ./  (  (f  0  . /f  req-f  req/f  0)  ./v2  +  l/z"2)  ) -amplit .  "2; 
rmserr=sum (err . A2)  ; 

%  FUNCTION  "ERR03.M"  FOR  "FMIN" .  LEAST  SQUARS  TECHNIQUE 
function  fplus=erro3 (wl, itrgh) 
[z, dens, s speed, vise] =argo3f (wl,  itrgh) ; 
gamma=5/3* (1+eps) ; 
npr=2/3* (1+eps)  ; 
sqri= ( (1+ j) /sqrt (2) ) * (1+eps) ; 
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f ofw=sqri*sqrt  (densA'3*sspeed'>4*npr .  .  . 
./(wl.*visc))./ ( gamma- 1) ; 

sl=real (fofw) ; 

s2=real  (z) ; 

s3=imag (fofw) ; 

s4=imag (z) ; 
fplus=(slA2-s2*2) "2+ (s3"2-s4~2) "2; 


%  FUNCTION  "ARG03F.M" .  "USED  BY  ERR03.M" 

function 
[z, dens, s speed, vise] =argo3 (wl, itrgh) 
w=wl; 
%   establish  some  often  used  constants 

s=5;  %  number  of  sections 

j=sqrt (-1) ; 
npr=2/3* (1+eps) ; 
gamma=5/3* (1+eps) ; 
cp= (2.5*8.314  3/ (4.e-3) ) * (1+eps) ; 
ksolid=. 16* (1+eps) ; 
sqri=( (1+j) /sqrt (2) ) * (1+eps) ; 
teren=zeros  (s, 1) ; 
numblay=zeros  (1, s) ; 
lengsec=zeros  (1, s) ; 
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tempR=zeros (1,  s)  ; 
tempL=zeros (1 ,  s) ; 
poros=zeros (1,  s)  ; 
ratio=zeros (1,  s) ; 
%  SET  VALUES 

termin=sum (abs (' free'  )  )  ; 
ampres= (1 . Ole+5) * (1+eps)  ; 
dramp= (1 .e-8) ; 


%  section  "1" 

terenl  =  ' opentu'  ; 

teren (1,1) =  sum (abs (terenl) )  ; 

numb lay (1, 1) =1; 

lengsec (1, 1)  =  (65. 34e-2) * (1+eps)  ; 

ternpRd,  1)  =293*  (1+eps)  ; 

tempL (1, 1) =293* (1+eps) ; 

ratio  (1, l)  =  (1.917e-2) * (1+eps)  ; 

poros  (1,1) =1* (1+eps)  ; 

ares (1,  : ) =pi* ratio (1,  :  ) *  (1+eps) ; 
%  section  "2" 

teren2  =  ' hexch'  ; 

teren  (1, 2) =sum (abs (teren2)  )  ; 

numblay (1, 2) =1; 

lengsec (1, 2)  =  ( . 82e-2) * (1+eps)  ; 

tempR(l,2)=2  93* (1+eps) ; 
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tempL(l,2)=293* (1+eps) ; 

ratio(l,2)=(5.08e-4)  *  (1+eps)  ; 

poros (1, 2) =. 666* (1+eps) / 
%  section  "3" 

teren3=' stack' ; 

teren (1,3) =sum (abs (teren3) ) ; 

numblay (1, 3) =11; 

lengsecd,  3)  =  (2.59e-2)  *  (1+eps)  ; 

tempRd,  3)  =itrgh*  (1+eps)  ; 

tempL(l,3)=2  93* (1+eps) ; 

ratio  (1, 3)=(8.636e-4) * (1+eps) ; 

poros (1,3) =.8  9473* (1+eps) ; 
%  section  "4" 

teren4  =  ' hexch'  ; 

teren (1, 4) =sum (abs (teren4) ) ; 

numblay (1, 4) =1; 

lengsec (1,4)= (.82e-2) * (1+eps) ; 

tempRd,  4)  =itrgh*  (1+eps)  ; 

tempL (1,4) =itrgh* (1+eps) ; 

ratio(l,4)=(5.08e-4) * (1+eps) ; 

poros (1, 4) =. 666* (1+eps) ; 
%  section  "5" 

teren5=' opentu' ; 

teren (1,5) =sum (abs (teren5) ) ; 

numblay (1, 5) =1; 

lengsec (1,5)= (. 6852) * (1+eps) ; 
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tempRd 
tempL (1 
ratio  (1 
poros  (1 
if  s>= 

%  section 
teren6= 
teren (1 
numblay 
lengsec 
tempR (1 
tempL (1 
ratio (1 
poros  (1 

%  section 
teren7= 
teren (1 
numblay 
lengsec 
tempRd 
tempL (1 
ratio  (1 
poros  (1 
end 


5) =itrgh* (1+eps) ; 
5) =itrgh* (1+eps) ; 
5)=(1. 917e-2) * (1+eps) ; 
5)=1* (1+eps) ; 

6" 

hexch' ; 

6) =sum (abs (teren6) ) ; 

1,6)=1; 

1, 6)=(1.02e-2) * (1+eps) ; 

6) =293*  (1+eps) ; 

6) =293*  (1+eps) ; 

6)=(1.02e-3) * (1+eps) ; 

6)=. 667* (1+eps) ; 


7" 


opentu' ; 

7) =sum (abs (teren7) ) ; 

1,7)=1; 

l,7)=(5.483e-2) * (1+eps) ; 

7)=293* (1+eps) ; 

7) =293*  (1+eps) ; 

7)=1.917e-2* (1+eps) ; 

7)=1* (1+eps); 


%  START  AT  THE  RIGHT  AND  MOVE  AT  THE  LEFT 
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numtot=sum (numblay) +1; 
dens (1, numtot) =ampres*4 . Oe-3  / (tempR(lf s) *8 . 3143) ;  %  density 
vised,  numtot)  =1.8  87e-5*  (tempR(l,s)  /273.15)~(.6567); 

%  viscosity 
kgas  (1, numtot ) =visc (1, numtot ) *cp/npr;  %  termal   conductivity 
sspeed (1, numtot) =972. 8*sqrt (tempR(l,s) /273. 15) ;  %   speed 

%   isothermal 
typel=' free' ; 
type2=' rigid' ; 

if  termin==sum (abs (type2) ) 
%   IMPEDANCE  OF  RIGID  TERMINATION, 
f ac3 ( 1 , numtot ) =sqrt ( (dens ( 1 , numtot ) . . . 

* sspeed (1, numtot) . A2) . / (w. *visc (1, numtot) ) ) ; 
z (1, numtot) = (1+ j) . * (dens (1, numtot) . *sspeed (1, numtot) . . . 

. *fac3 (1, numtot) *sqrt (npr) ) / (sqrt (2) * ( gamma- 1) ) ; 
%  IMPEDANCE  OF  FREE  TERMINATION 
else 
z (1, numtot) =f ree (dens (1, numtot) , vise (1, numtot) , w, ratio (1, s) . 
, sspeed (1, numtot) , gamma, npr) ; 
end 

%  IMPEDANCE  TRANSLATE  FOR  THE  OPEN  TUBE  OR  HEAT  EXCHANGER. 
fault=' stack' ; 
truel=' opentu' ; 
true2=' hexch' ; 
for  k=s:-l:l 
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jup=0; 

jup=numtot-sum ( jup+numblay (1 , k : s) ) / 
if  teren (1, k) ~=sum (abs (fault ) ) 
dens (1, jup) =ampres*4 . Oe-3  / (tempR (1, k) . *8 . 3143) ;  %  density 
vise (1, jup)=1.887e-5.* (tempR(l,k) /273.15)  .  A  (.65  67) ; 

%  viscosity 
kgas (1, jup) =visc (1, jup) . *cp/npr;   %  termal   conductivity 
sspeedd,  jup) =972 . 8 . *sqrt  (tempR(l,k)  /273.15)  ;    %   speed 

%   isothermal 

%      GET  LAMBDA  AND  LAMBDA (T) 
lambda (1, jup) =ratio  (1, k)  . *sqrt (dens (1, jup)  . *w. /vise (1, jup) ) ; 
lambdt (1, jup) =sqrt (npr) . *lambda (1, jup) ; 
%      GET  F(L)/  F(L(T))  AND  WAVENUMBERS  FOR  OPEN  TUBE 
PARTS 

if  teren (1, k) ==sum (abs (truel) ) 
flam(l, jup)=l-(l+j) .*sqrt (2) . /lambda (1, jup) ;    %       f(l) 
flamt (1, jup)=l-(l+j) .*sqrt (2) ./lambdt (1, jup) ;   %  f(l(t)) 
fac31= (1+ ( gamma- 1) /sqrt (npr) ) /sqrt (2) ; 
ka (1, jup) = (w. /s speed (1, jup) ) . . . 

. *(l+(l+j) . *fac31. /lambda (1, jup) ) ; 
else 
%      GET  F(L),  F(L(T))  AND  WAVENUMBERS  FOR  HEAT 
EXCHANGERS  TYPE  "SLIT". 

sqrmi= (1- j) /sqrt (2) ; 

ar (1, jup) =sqrmi . *  lambda (1, jup) ; 
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argum (1, jup) =exp (-ar (1, jup) ) ; 

ar (1, jup) =ar (1, jup) /2; 

ctanh (1, jup) = (1-argum (1, jup) ) . / (1+argum (1, jup) ) ; 

flam(l, jup) =l-ctanh (1, jup) . /ar (1, jup) ;    %  f(l) 

ar (1, jup) =sqrmi*lambdt (1, jup) ; 

argumd,  jup)  =exp  (-ar  (1,  jup)  )  ; 

ar  (1, jup)=ar (1, jup) 12; 

ctanh (1, jup) = (1-argum (1, jup) ) . / (1+argum (1, jup) ) ; 

flamt (1,  jup)=l-ctanh(l, jup) ./ar (1, jup) ;  %  f(l(t)) 

low=w/sspeed ( 1 , jup) ; 

ka (1, jup) =low. *sqrt ( (gamma- ( gamma- 1) . . . 

. *flamt (1,  jup) ) . /flam(l, jup) ) ; 
end 
%  TRANSLATION  THEOREM 
comden (1, jup) =dens (1, jup) . /flam (1, jup) ; 
zint (1, jup) =comden (1, jup) . *w. / (ka (1, jup) . *poros (1, k) ) ; 
dsub (1, k) =lengsec  (1, k)  . /numblay (1, k) ; 
ar (1,  jup) =ka (1, jup)  . *dsub (1, k) ; 
sn (1,  jup) =sin (ar (1, jup) ) ; 
cs (1,  jup) =cos (ar (1, jup) ) ; 
ct (1,  jup) =cs (1, jup)  . /sn (1, jup) ; 
fac3 (1,  jup) =zint (1, jup) . /z (1, 1+ jup) ; 
z (1,  jup) =zint (1, jup)  . * (ct (1, jup)  . . . 

-j.*fac3(l, jup) ) ./ (fac3(l, jup) . *ct (1, jup)-j) ; 
else 
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%   STACK  WITH  LINEAR  DISTRIBUTION  BETWEEN  THE  ENDS  OF  THE 
STACK. 

toz= (tempR(l, k) -tempL (1, k) ) /lengsec (1, k) ;    %  approx 
dz/dt 

ns=0; 

for  la=k+numblay (1, k) -1 : -1 : k 

ns=ns+l; 

tens=tempR(l, k) - (tempR (1, k) -tempL (1,  k) ) *ns/numblay (1,  k) ; 
tnsml=tempR(l, k) - (tempR(l, k) -tempL (1, k) ) * (ns-1) /numblay (1, k) 

tave (1, la) = (tens+tnsml) 12; 

dsub (1, k) = lengsec (1, k) /numblay (1, k) ; 

dens  (1, la)=ampres*4.0e-3  /  (tave (1, la) *8 . 3143) ;     %  density 

visc(l,la)=1.887e-5* (tave(l, la) /273. 15) A ( . 6567) ;   % 
viscosity 

kgas (1, la) =visc (1, jup) *cp/npr;  %  termal  conductivity 
sspeed(l,la)=972.8*sqrt (tave(l,la)/273.15) ;       %   speed 

%   isothermal 
%      GET  LAMBDA  AND  LAMBDA (T) 
lambda (1, la) =ratio (1, k) . *sqrt (dens (1,1a) .*w./visc(l,la)); 

lambdt (1, la) =sqrt (npr) *  lambda (1, la) ; 

ar  (1, la) =sqrmi* lambda (1,1a); 

argum (l,la)=exp(-ar(l,la)  ); 

ar (1, la)=ar (1, la) /2; 

ctanh (1,  la)  =  (1 -argum (1, la) )  . / (1  +  argum (1, la) ) ; 

flam(l,la)=l-ctanh(l, la)  . /ar  (1, la) ;        %  f(l) 
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ar (1,  la) =sqrmi*lambdt (1, la) ; 

argum (l,la)=exp(-ar(l,la) ) ; 

ar  (l,la)=ar (1,1a) 12; 

ctanh  (1, la)  =  (1 -argum (1, la) )  . / (1+argum (1,  la)  )  ; 

flamt (1, la)=l-ctanh(l, la)  . /ar  (1, la) ;       %  f(l(t)) 

low=w/sspeed (1, la) ; 

ka (1, la) =low. *sqrt ( (gamma- (gamma-1) . . . 

. * flamt (1, la) ) . /flam(l, la) ) ; 
zint (1, la) =dens (1,1a) .  *w. . . 

./ (poros (l,k) .*f lam (1,1a) . *ka(l,la) ) ; 
alprimd,  la) =toz*  (flamt  (1,  la)  .  .  . 

. /flam (1,1a) -1)  ./  (2*tave (1, la) * (1-npr) ) ; 
kal=ka(l,la) . ' ; 
zintl  =  zint (1, la) . '  ; 
alprl=alprim (1, la) . '; 

zetaO= (lengsec  (1, k) - (ns-1) *lengsec (1, k) /numblay (1, k) ) ; 
zetaf in=  (lengsec  (1, k) - (ns) * lengsec (1, k) /numblay (1, k) ) ; 
zO=z (1, 1+la) .' ; 
tol=l.e-4; 
[zeta, zout] =t rode 4 5 (' tube' ,zetaO,zetafin,zO,tol,kal,zintl,al 
prl); 

n=length (zout) ; 
z (1,  la) =zout (n, 1)  . ' ; 
end 
end 
end 
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2  =  2(1,1)  ; 
dens=dens (1,1); 
sspeed=sspeed (1,1); 
visc=visc  (1,1); 
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